/* Finite element subroutines: mesh data structures
   Jukka Toivanen, Jan Stebel
*/

#include <math.h>
#include <map>
#include <list>
#include <algorithm>
#include <cstdarg>
#include <string.h>
#include "mesh.h"



void intersection(ChildParentListType &set1, ChildParentListType &set2, ChildParentListType &result)
{
  ChildParentListType::iterator first1 = set1.begin(),
                      last1  = set1.end(),
                      first2 = set2.begin(),
                      last2  = set2.end();
  while (first1!=last1 && first2!=last2)
  {
    if (*first1<*first2) ++first1;
    else if (*first2<*first1) ++first2;
    else { result.push_back(*first1++); first2++; }
  }
}


CElement::CElement(const CElement &e)
    : nphysical(e.nphysical),
      ngeom(e.ngeom),
      npart(e.npart),
      type(e.type)
{
  children = e.children;
  parents = e.parents;
  parentNodes = e.parentNodes;
  nodes = e.nodes;
}



int Mesh::addMidNode(int size, ...)
{
  CVar x = 0, y = 0, z = 0;
  va_list node;
  
  va_start(node, size);
  for (int i=0; i<size; i++)
  {
    int n = va_arg(node, int);
    x += coords[3*n  ];
    y += coords[3*n+1];
    z += coords[3*n+2];
  }
  va_end(node);
  coords.push_back(x/size);
  coords.push_back(y/size);
  coords.push_back(z/size);
  nnodes++;
  return nnodes-1;
}


void Mesh::splitVolume(int i)
{
  CElement tet(volumes[i]);
  vector<int> node;
  int new_node;
  
  for (unsigned int j=0; j<volumes[i].nodes.size(); j++)
    node.push_back(volumes[i].nodes[j]);
  
  // add midpoint
  new_node = addMidNode(4, node[0], node[1], node[2], node[3]);
  
  switch (volumes[i].type)
  {
  case GMSH_TETRAHEDRON:
      
      // 0-1-2-4
      volumes[i].nodes[3] = new_node;
      // 0-1-4-3
      tet.nodes[2] = new_node;
      tet.nodes[3] = node[3];
      volumes.push_back(tet);
      // 0-2-3-4
      tet.nodes[1] = node[2];
      tet.nodes[2] = node[3];
      tet.nodes[3] = new_node;
      volumes.push_back(tet);
      // 1-2-4-3
      tet.nodes[0] = node[1];
      tet.nodes[2] = new_node;
      tet.nodes[3] = node[3];
      volumes.push_back(tet);
      break;
  case GMSH_QUAD_TETRAHEDRON:
      // add point between 0 and midpoint
      int node0n = addMidNode(2, node[0], new_node);
      // add point between 1 and midpoint
      int node1n = addMidNode(2, node[1], new_node);
      // add point between 2 and midpoint
      int node2n = addMidNode(2, node[2], new_node);
      // add point between 3 and midpoint
      int node3n = addMidNode(node[3], new_node);
      // update linear nodes
      lnmap.resize(nnodes);
      lnmap[new_node] = nlnodes;
      lnmap[node0n] = -1;
      lnmap[node1n] = -1;
      lnmap[node2n] = -1;
      lnmap[node3n] = -1;
      nlnodes++;
      
      // 0-1-2-N
      volumes[i].nodes[3] = new_node;
      volumes[i].nodes[7] = node0n;
      volumes[i].nodes[8] = node2n;
      volumes[i].nodes[9] = node1n;
      // 0-1-N-3
      tet.nodes[0] = node[0];
      tet.nodes[1] = node[1];
      tet.nodes[2] = new_node;
      tet.nodes[3] = node[3];
      tet.nodes[4] = node[4];
      tet.nodes[5] = node1n;
      tet.nodes[6] = node0n;
      tet.nodes[7] = node[7];
      tet.nodes[8] = node3n;
      tet.nodes[9] = node[9];
      volumes.push_back(tet);
      // 0-2-3-N
      tet.nodes[0] = node[0];
      tet.nodes[1] = node[2];
      tet.nodes[2] = node[3];
      tet.nodes[3] = new_node;
      tet.nodes[4] = node[6];
      tet.nodes[5] = node[8];
      tet.nodes[6] = node[7];
      tet.nodes[7] = node0n;
      tet.nodes[8] = node3n;
      tet.nodes[9] = node2n;
      volumes.push_back(tet);
      // 1-2-N-3
      tet.nodes[0] = node[1];
      tet.nodes[1] = node[2];
      tet.nodes[2] = new_node;
      tet.nodes[3] = node[3];
      tet.nodes[4] = node[5];
      tet.nodes[5] = node2n;
      tet.nodes[6] = node1n;
      tet.nodes[7] = node[9];
      tet.nodes[8] = node3n;
      tet.nodes[9] = node[8];
      volumes.push_back(tet);
      break;
  }
}


void Mesh::splitFace(int i)
{
  CElement e(faces[i]);
  vector<int> node;
  int new_node;
  
  for (unsigned int j=0; j<faces[i].nodes.size(); j++)
    node.push_back(faces[i].nodes[j]);
  
  // add midpoint
  new_node = addMidNode(3, node[0], node[1], node[2]);
  
  switch (faces[i].type)
  {
  case GMSH_TRIANGLE:
      // 0-1-N
      faces[i].nodes[2] = new_node;
      // 0-N-2
      e.nodes[1] = new_node;
      faces.push_back(e);
      // N-1-2
      e.nodes[0] = new_node;
      e.nodes[1] = node[1];
      faces.push_back(e);
      break;
  case GMSH_QUAD_TRIANGLE:
      // add point between 0 and midpoint
      int node0n = addMidNode(2, node[0], new_node);
      // add point between 1 and midpoint
      int node1n = addMidNode(2, node[1], new_node);
      // add point between 2 and midpoint
      int node2n = addMidNode(2, node[2], new_node);
      // update linear nodes
      lnmap.resize(nnodes);
      lnmap[new_node] = nlnodes;
      lnmap[node0n] = -1;
      lnmap[node1n] = -1;
      lnmap[node2n] = -1;
      nlnodes++;
      
      // 0-1-N
      faces[i].nodes[2] = new_node;
      faces[i].nodes[4] = node1n;
      faces[i].nodes[5] = node0n;
      // 0-N-2
      e.nodes[1] = new_node;
      e.nodes[3] = node0n;
      e.nodes[4] = node2n;
      faces.push_back(e);
      // N-1-2
      e.nodes[0] = new_node;
      e.nodes[1] = node[1];
      e.nodes[3] = node1n;
      e.nodes[4] = node[4];
      e.nodes[5] = node2n;
      faces.push_back(e);
      break;
  }
}


void Mesh::destroy()
{
  coords.clear();
  physicalNames.clear();
  volumes.clear();
  faces.clear();
  edges.clear();
  points.clear();
  lnmap.clear();
}

Mesh::~Mesh()
{
  destroy();
}

void Mesh::deform(CVar& func(double, double, double, int))
{
  int i;
  double x, y, z;
  
  // loop through all nodes
  for (i=0; i<nnodes; i++)
  {
    x = coords[3*i  ].getVal();
    y = coords[3*i+1].getVal();
    z = coords[3*i+2].getVal();
    
    coords[3*i  ] = func(x, y, z, 0);
    coords[3*i+1] = func(x, y, z, 1);
    coords[3*i+2] = func(x, y, z, 2);
  }
}

CVar& Mesh::getCoord(int node, int coord)
{
  return coords[3*node+coord];
}

void Mesh::setCoord(int node, int coord, CVar x)
{
  coords[3*node+coord] = x;
}

void Mesh::calculateFaceNormal(int fid, CVec &n, CVec &a, CVec &b)
{
  int i;
  double anorm, bnorm, nnorm;
  a.resize(3);
  b.resize(3);
  n.resize(3);
  for (i=0; i<3; i++)
  {
    a(i) = getCoord(getFaceNode(fid, 1), i) - getCoord(getFaceNode(fid, 0), i);
    b(i) = getCoord(getFaceNode(fid, 2), i) - getCoord(getFaceNode(fid, 0), i);
  }
  anorm = a.norm2();
  bnorm = b.norm2();
  for (i=0; i<3; i++)
  {
    a(i) = a(i)/anorm;
    b(i) = b(i)/bnorm;
  }
  n(0) = a(1)*b(2) - a(2)*b(1);
  n(1) = a(2)*b(0) - a(0)*b(2);
  n(2) = a(0)*b(1) - a(1)*b(0);
  nnorm = n.norm2();
  for (i=0; i<3; i++) n(i) = n(i)/nnorm;
}

void Mesh::calculateEdgeTangent(int eid, CVec &t)
{
  int i;
  double norm = 0;
  t.resize(3);
  for (i=0; i<3; i++)
  {
    t(i) = getCoord(getEdgeNode(eid, 1), i) - getCoord(getEdgeNode(eid, 0), i);
    norm += t(i).getVal()*t(i).getVal();
  }
  norm = sqrt(norm);
  for (i=0; i<3; i++) t(i) = t(i)/norm;
}



void Mesh::readGMSH(string ifname, int partition)
{
  int dim, i, j, ip, nelems, n, ntags, nparts, tag;
  double x, y, z;
  char version[255];
  string s, s1;
  ifstream f(ifname.c_str());
  map<int, int> mapNodes;
  CElement el;
  
  fname = ifname;

  try
  {
    getline(f, s);
    if (s.compare("$MeshFormat") != 0) throw 1;
    getline(f, s);
    sscanf(s.c_str(), "%s 0 %d", version, &i);
    if (i != sizeof(double)) throw 2;
    getline(f, s);
    if (s.compare("$EndMeshFormat") != 0) throw 1;
  
    getline(f, s);
    if (s.compare("$PhysicalNames") == 0)
    { // read through the physical names block
      getline(f, s);
      sscanf(s.c_str(), "%d", &n);
      
      for (i=0; i<n; i++)
      {
        getline(f, s);
        if (strcmp(version, "2.0") == 0 || strcmp(version, "2") == 0)
        {
          sscanf(s.c_str(), "%d", &ip);
        }
        else
        {
          sscanf(s.c_str(), "%d %d", &dim, &ip);
        }
	j  = s.find("\"");
	s1 = s.substr(j+1, s.length()-j-2);
	physicalNames[ip] = s1;
	physNameIds[s1] = ip;
      }
      getline(f, s); // end tag "$EndPhysicalNames"
      getline(f, s);
    }
    if (s.compare("$Nodes") == 0)
    { // read nodes data
      getline(f, s);
      sscanf(s.c_str(), "%d", &nnodes);
      coords.resize(3*nnodes);
      lnmap.resize(nnodes);
      node_partition.resize(nnodes);
      nlnodes = 0;
      for (i=0; i<nnodes; lnmap[i] = -1, i++);

      for (i=0; i<nnodes; i++)
      {
        getline(f, s);
        sscanf(s.c_str(), "%d %lf %lf %lf", &j, &x, &y, &z);
        mapNodes[j] = i;
        setCoord(i, 0, x);
        setCoord(i, 1, y);
        setCoord(i, 2, z);
        node_partition[i] = 0;
      }
  
      getline(f, s);
      sscanf(s.c_str(), "$EndNodes");
    }
    else
    {
      throw 3;
    }
  
    getline(f, s);
    if (s.compare("$Elements") != 0) throw 3;

    getline(f, s);
    sscanf(s.c_str(), "%d", &nelems);
  
    for (int i=0; i<nelems; i++)
    {
      bool shared_element = false,
    	   element_of_this_partition = false;
      f >> n >> el.type >> ntags;
      if (ntags >= 3) // read tags
      {
        f >> el.nphysical;
        f >> el.ngeom;
        f >> nparts;
        if (ntags >= 4)
        {
        	f >> el.npart;

        	if (nparts > 1) shared_element = true;

        	for (j=1; j<nparts; j++)
        	{
        		f >> tag;
        		if (fabs(tag) == partition) element_of_this_partition = true;
        	}

        	for (j=3+nparts; j<ntags; f >> tag, j++);
        }
      }
      else if (ntags >= 2) // read tags
      {
        f >> el.nphysical;
        f >> el.ngeom;
        el.npart = 0;
        for (j=2; j<ntags; f >> tag, j++);
      }
      else if (ntags >= 1) // read tags
      {
        f >> el.nphysical;
        el.ngeom = 0;
        el.npart = 0;
        for (j=1; j<ntags; f >> tag, j++);
      }
      else
      {
        el.npart     = 0;
        el.ngeom     = 0;
        el.nphysical = 0;
        for (j=0; j<ntags; f >> tag, j++);
      }

      el.nodes.clear();
      int el_nodes=0, el_lnodes=0;
      switch (el.type)
      {
        case GMSH_TETRAHEDRON:
          el_nodes  = 4;
          el_lnodes = 4;
          break;
        case GMSH_TRIANGLE:
          el_nodes  = 3;
          el_lnodes = 3;
          break;
        case GMSH_LINE:
          el_nodes  = 2;
          el_lnodes = 2;
          break;
        case GMSH_POINT:
          el_nodes  = 1;
          el_lnodes = 1;
	  break;
        case GMSH_QUAD_TETRAHEDRON:
          el_nodes  = 10;
          el_lnodes = 4;
          break;
        case GMSH_QUAD_TRIANGLE:
          el_nodes  = 6;
          el_lnodes = 3;
          break;
        case GMSH_QUAD_LINE:
          el_nodes  = 3;
          el_lnodes = 2;
          break;
      }

      if (partition_nodes.find(el.npart) == partition_nodes.end())
      {
    	  set<int> v;
    	  v.clear();
    	  partition_nodes[el.npart] = v;
      }

      for (j=0; j<el_nodes; j++) // read nodes
      {
        f >> n;
        el.nodes.push_back(mapNodes[n]);

        // mark nodes that belong to more partitions
        if (shared_element || (partition != 0 && node_partition[mapNodes[n]] != 0 && node_partition[mapNodes[n]] != el.npart))
        	shared_nodes.insert(mapNodes[n]);

        if (partition == 0)
        {
        	partition_nodes[0].insert(mapNodes[n]);
        }
        else if (node_partition[mapNodes[n]] == 0 || el.npart < node_partition[mapNodes[n]])
        {
        	partition_nodes[node_partition[mapNodes[n]]].erase(mapNodes[n]);
        	node_partition[mapNodes[n]] = el.npart;
        	partition_nodes[el.npart].insert(mapNodes[n]);
        }
      }
      for (j=0; j<el_lnodes; j++) // mark linear nodes
      {
        if (lnmap[el.nodes[j]]==-1)
        {
          lnmap[el.nodes[j]] = nlnodes;
          nlnodes++;
        }
      }

      if (element_of_this_partition)
    	  el.npart = partition;

      // skip elements that do not belong to the given partition
      if (partition != 0 && el.npart != partition) continue;

      switch (el.type)
      {
      case GMSH_TETRAHEDRON:
      case GMSH_QUAD_TETRAHEDRON:
          volumes.push_back(el);
          break;
      case GMSH_TRIANGLE:
      case GMSH_QUAD_TRIANGLE:
          faces.push_back(el);
          break;
      case GMSH_LINE:
      case GMSH_QUAD_LINE:
          edges.push_back(el);
          break;
      case GMSH_POINT:
          points.push_back(el);
          break;
      }
    }
  }
  catch (int e)
  {
    switch (e)
    {
      case 1: cout << "Unsupported mesh file format.\n"; break;
      case 2: cout << "Data size does not match this system.\n"; break;
      case 3: cout << "Error reading mesh file.\n"; break;
      case 4: cout << "Element can belong to only one mesh partition."; break;
    }
    f.close();
    mapNodes.clear();
    exit(1);
  }
  f.close();
  mapNodes.clear();
}


void Mesh::splitBoundaryVolumes()
{
  bool isBoundaryNode[nnodes];
  bool boundaryVolume;
  int i, j;
  
  for (i=0; i<nnodes; i++) isBoundaryNode[i] = false;
  
  // mark boundary nodes
  for (i=0; i<(int)faces.size(); i++)
    for (j=0; j<(int)faces[i].nodes.size(); j++)
      isBoundaryNode[faces[i].nodes[j]] = true;

  // find volumes whose all vertices lie at the boundary
  int nvolumes = volumes.size();
  for (i=0; i<nvolumes; i++)
  {
    boundaryVolume = true;
    for (j=0; j<(int)volumes[i].nodes.size(); j++)
    {
      if (!isBoundaryNode[volumes[i].nodes[j]])
      {
        boundaryVolume = false;
	break;
      }
    }
    // split the volume
    if (boundaryVolume) splitVolume(i);
  }
}







int LinMesh::addEdge(int n1, int n2, int nphys, int ngeom, int npart)
{
  CElement e;
  e.nodes.push_back(n1);
  e.nodes.push_back(n2);
  e.nphysical = nphys;
  e.ngeom = ngeom;
  e.npart = npart;
  edges.push_back(e);
  return (int)edges.size()-1;
}

int LinMesh::addFace(int n1, int n2, int n3, int nphys, int ngeom, int npart)
{
  CElement e;
  e.nodes.push_back(n1);
  e.nodes.push_back(n2);
  e.nodes.push_back(n3);
  e.nphysical = nphys;
  e.ngeom = ngeom;
  e.npart = npart;
  faces.push_back(e);
  return (int)faces.size()-1;
}







void LinMesh::findChildrenParents(bool create_new)
{
  map<int,ChildParentListType> nodeFaces;
  map<int,ChildParentListType> nodeEdges;

  // each entry in nodeFaces array contains a list of faces that use the particular node
  for (unsigned int i=0; i<faces.size(); i++)
  {
    for (int j=0; j<getNfaceDofs(); j++)
    {
      nodeFaces[faces[i].nodes[j]].push_back(i);
    }
  }
  // the same for edges
  for (unsigned int i=0; i<edges.size(); i++)
  {
    for (int j=0; j<getNedgeDofs(); j++)
    {
      nodeEdges[edges[i].nodes[j]].push_back(i);
    }
  }
  
  // sort lists
  for (int i=0; i<nnodes; i++)
  {
    sort(nodeFaces[i].begin(), nodeFaces[i].end());
    sort(nodeEdges[i].begin(), nodeEdges[i].end());
  }
  
  // for each volume find its faces
  for (unsigned int i=0; i<volumes.size(); i++)
  {
    ChildParentListType result;
    map<int, int> nodeVolumeOrder;
    
    for (int j=0; j<getNvolumeDofs(); j++) nodeVolumeOrder[getVolumeNode(i,j)] = j;
    
    volumes[i].children.clear();
    
    result.clear();
    intersection(nodeFaces[volumes[i].nodes[0]], nodeFaces[volumes[i].nodes[1]], result);
    intersection(result, nodeFaces[volumes[i].nodes[2]], volumes[i].children);
    if (create_new && volumes[i].children.size() != 1)
    {
      int nf = addFace(volumes[i].nodes[0], volumes[i].nodes[1], volumes[i].nodes[2]);
      nodeFaces[volumes[i].nodes[0]].push_back(nf);
      nodeFaces[volumes[i].nodes[1]].push_back(nf);
      nodeFaces[volumes[i].nodes[2]].push_back(nf);
      sort(nodeFaces[volumes[i].nodes[0]].begin(), nodeFaces[volumes[i].nodes[0]].end());
      sort(nodeFaces[volumes[i].nodes[1]].begin(), nodeFaces[volumes[i].nodes[1]].end());
      sort(nodeFaces[volumes[i].nodes[2]].begin(), nodeFaces[volumes[i].nodes[2]].end());
      volumes[i].children.push_back(nf);
    }
    
    result.clear();
    intersection(nodeFaces[volumes[i].nodes[0]], nodeFaces[volumes[i].nodes[1]], result);
    intersection(result, nodeFaces[volumes[i].nodes[3]], volumes[i].children);
    if (create_new && volumes[i].children.size() != 2)
    {
      int nf = addFace(volumes[i].nodes[0], volumes[i].nodes[1], volumes[i].nodes[3]);
      nodeFaces[volumes[i].nodes[0]].push_back(nf);
      nodeFaces[volumes[i].nodes[1]].push_back(nf);
      nodeFaces[volumes[i].nodes[3]].push_back(nf);
      sort(nodeFaces[volumes[i].nodes[0]].begin(), nodeFaces[volumes[i].nodes[0]].end());
      sort(nodeFaces[volumes[i].nodes[1]].begin(), nodeFaces[volumes[i].nodes[1]].end());
      sort(nodeFaces[volumes[i].nodes[3]].begin(), nodeFaces[volumes[i].nodes[3]].end());
      volumes[i].children.push_back(nf);
    }
    
    result.clear();
    intersection(nodeFaces[volumes[i].nodes[0]], nodeFaces[volumes[i].nodes[2]], result);
    intersection(result, nodeFaces[volumes[i].nodes[3]], volumes[i].children);
    if (create_new && volumes[i].children.size() != 3)
    {
      int nf = addFace(volumes[i].nodes[0], volumes[i].nodes[2], volumes[i].nodes[3]);
      nodeFaces[volumes[i].nodes[0]].push_back(nf);
      nodeFaces[volumes[i].nodes[2]].push_back(nf);
      nodeFaces[volumes[i].nodes[3]].push_back(nf);
      sort(nodeFaces[volumes[i].nodes[0]].begin(), nodeFaces[volumes[i].nodes[0]].end());
      sort(nodeFaces[volumes[i].nodes[2]].begin(), nodeFaces[volumes[i].nodes[2]].end());
      sort(nodeFaces[volumes[i].nodes[3]].begin(), nodeFaces[volumes[i].nodes[3]].end());
      volumes[i].children.push_back(nf);
    }
    
    result.clear();
    intersection(nodeFaces[volumes[i].nodes[1]], nodeFaces[volumes[i].nodes[2]], result);
    intersection(result, nodeFaces[volumes[i].nodes[3]], volumes[i].children);
    if (create_new && volumes[i].children.size() != 4)
    {
      int nf = addFace(volumes[i].nodes[1], volumes[i].nodes[2], volumes[i].nodes[3]);
      nodeFaces[volumes[i].nodes[1]].push_back(nf);
      nodeFaces[volumes[i].nodes[2]].push_back(nf);
      nodeFaces[volumes[i].nodes[3]].push_back(nf);
      sort(nodeFaces[volumes[i].nodes[1]].begin(), nodeFaces[volumes[i].nodes[1]].end());
      sort(nodeFaces[volumes[i].nodes[2]].begin(), nodeFaces[volumes[i].nodes[2]].end());
      sort(nodeFaces[volumes[i].nodes[3]].begin(), nodeFaces[volumes[i].nodes[3]].end());
      volumes[i].children.push_back(nf);
    }
    
    for (ChildParentListType::iterator it=volumes[i].children.begin(); it!=volumes[i].children.end(); it++)
    {
      faces[*it].parents.push_back(i);
      
      for (int j=0; j<getNfaceDofs(); j++)
      {
        faces[*it].parentNodes[j][faces[*it].parents.size()-1] = nodeVolumeOrder[getFaceNode(*it,j)];
      }
    }
    nodeVolumeOrder.clear();
  }
  // for each face find its edges
  for (unsigned int i=0; i<faces.size(); i++)
  {
    map<int, int> nodeFaceOrder;
    
    faces[i].children.clear();
    intersection(nodeEdges[faces[i].nodes[0]], nodeEdges[faces[i].nodes[1]], faces[i].children);
    if (create_new && faces[i].children.size() != 1)
    {
      int ne = addEdge(faces[i].nodes[0], faces[i].nodes[1]);
      nodeEdges[faces[i].nodes[0]].push_back(ne);
      nodeEdges[faces[i].nodes[1]].push_back(ne);
      sort(nodeEdges[faces[i].nodes[0]].begin(), nodeEdges[faces[i].nodes[0]].end());
      sort(nodeEdges[faces[i].nodes[1]].begin(), nodeEdges[faces[i].nodes[1]].end());
      faces[i].children.push_back(ne);
    }
    intersection(nodeEdges[faces[i].nodes[0]], nodeEdges[faces[i].nodes[2]], faces[i].children);
    if (create_new && faces[i].children.size() != 2)
    {
      int ne = addEdge(faces[i].nodes[0], faces[i].nodes[2]);
      nodeEdges[faces[i].nodes[0]].push_back(ne);
      nodeEdges[faces[i].nodes[2]].push_back(ne);
      sort(nodeEdges[faces[i].nodes[0]].begin(), nodeEdges[faces[i].nodes[0]].end());
      sort(nodeEdges[faces[i].nodes[2]].begin(), nodeEdges[faces[i].nodes[2]].end());
      faces[i].children.push_back(ne);
    }
    intersection(nodeEdges[faces[i].nodes[1]], nodeEdges[faces[i].nodes[2]], faces[i].children);
    if (create_new && faces[i].children.size() != 3)
    {
      int ne = addEdge(faces[i].nodes[1], faces[i].nodes[2]);
      nodeEdges[faces[i].nodes[1]].push_back(ne);
      nodeEdges[faces[i].nodes[2]].push_back(ne);
      sort(nodeEdges[faces[i].nodes[1]].begin(), nodeEdges[faces[i].nodes[1]].end());
      sort(nodeEdges[faces[i].nodes[2]].begin(), nodeEdges[faces[i].nodes[2]].end());
      faces[i].children.push_back(ne);
    }
    
    for (int j=0; j<getNfaceDofs(); j++) nodeFaceOrder[getFaceNode(i,j)] = j;
    
    for (ChildParentListType::iterator it=faces[i].children.begin(); it!=faces[i].children.end(); it++)
    {
      edges[*it].parents.push_back(i);
      
      for (int j=0; j<getNedgeDofs(); j++)
      {
        edges[*it].parentNodes[j][edges[*it].parents.size()-1] = nodeFaceOrder[getEdgeNode(*it,j)];
      }
    }
    nodeFaceOrder.clear();
  }
  
  nodeFaces.clear();
  nodeEdges.clear();
}









void QuadMesh::writeGMSH(string fname)
{
  ofstream f(fname.c_str());
  
  f << "$MeshFormat\n";
  f << "2 0 " << sizeof(double) << endl;
  f << "$EndMeshFormat\n";
  if (getNphysicalNames() > 0)
  {
    f << "$PhysicalNames\n";
    f << getNphysicalNames() << endl;
    for (int i=0; i<getNphysicalNames(); i++)
    {
      f << i+1 << " \"" << physicalNames[i] << "\"\n";
    }
    f << "$EndPhysicalNames\n";
  }
  f << "$Nodes\n";
  f << nnodes << endl;
  for (int i=0; i<nnodes; i++)
  {
    f << i+1 << " ";
    for (int j=0; j<3; j++)
    {
      f << getCoord(i, j).getVal();
      if (j<2) f << " ";
    }
    f << endl;
  }
  f << "$EndNodes\n";
  f << "$Elements\n";
  f << faces.size() + volumes.size() << endl;
  int e=0;
  for (unsigned int i=0; i<faces.size(); i++)
  {
    e++;
    f << e << " 9 3 " << getFacePhysicalNo(i) << " "
      << getFaceGeomNo(i) << " "
      << getFacePartNo(i) << " ";
    for (int j=0; j<6; j++)
    {
      f << getFaceNode(i, j)+1;
      if (j<5) f << " ";
    }
    f << endl;
  }
  for (unsigned int i=0; i<volumes.size(); i++)
  {
    e++;
    f << e << " 11 3 " << getVolumePhysicalNo(i) << " "
      << getVolumeGeomNo(i) << " "
      << getVolumePartNo(i) << " ";
    for (int j=0; j<10; j++)
    {
      f << getVolumeNode(i, j)+1;
      if (j<9) f << " ";
    }
    f << endl;
  }
  f << "$EndElements\n";
  f.close();
}



int QuadMesh::addEdge(int n1, int n2, int n3, int nphys, int ngeom, int npart)
{
  CElement e;
  e.nodes.push_back(n1);
  e.nodes.push_back(n2);
  e.nodes.push_back(n3);
  e.nphysical = nphys;
  e.ngeom = ngeom;
  e.npart = npart;
  edges.push_back(e);
  return (int)edges.size()-1;
}

void QuadMesh::findChildrenParents(bool create_new)
{
  map<int,ChildParentListType> nodeFaces;
  map<int,ChildParentListType> nodeEdges;

  // each entry in nodeFaces array contains a list of faces that use the particular node
  for (unsigned int i=0; i<faces.size(); i++)
  {
    for (int j=0; j<getNfaceDofs(); j++)
    {
      nodeFaces[faces[i].nodes[j]].push_back(i);
    }
  }
  // the same for edges
  for (unsigned int i=0; i<edges.size(); i++)
  {
    for (int j=0; j<getNedgeDofs(); j++)
    {
      nodeEdges[edges[i].nodes[j]].push_back(i);
    }
  }
  
  // sort lists
  for (int i=0; i<nnodes; i++)
  {
    sort(nodeFaces[i].begin(), nodeFaces[i].end());
    sort(nodeEdges[i].begin(), nodeEdges[i].end());
  }
  
  // for each volume find its faces
  for (unsigned int i=0; i<volumes.size(); i++)
  {
    ChildParentListType result;
    map<int, int> nodeVolumeOrder;
    
    for (int j=0; j<getNvolumeDofs(); j++) nodeVolumeOrder[getVolumeNode(i,j)] = j;
    
    volumes[i].children.clear();
    
    result.clear();
    intersection(nodeFaces[volumes[i].nodes[0]], nodeFaces[volumes[i].nodes[1]], result);
    intersection(result, nodeFaces[volumes[i].nodes[2]], volumes[i].children);
    
    result.clear();
    intersection(nodeFaces[volumes[i].nodes[0]], nodeFaces[volumes[i].nodes[1]], result);
    intersection(result, nodeFaces[volumes[i].nodes[3]], volumes[i].children);
    
    result.clear();
    intersection(nodeFaces[volumes[i].nodes[0]], nodeFaces[volumes[i].nodes[2]], result);
    intersection(result, nodeFaces[volumes[i].nodes[3]], volumes[i].children);
    
    result.clear();
    intersection(nodeFaces[volumes[i].nodes[1]], nodeFaces[volumes[i].nodes[2]], result);
    intersection(result, nodeFaces[volumes[i].nodes[3]], volumes[i].children);
    
    for (ChildParentListType::iterator it=volumes[i].children.begin(); it!=volumes[i].children.end(); it++)
    {
      faces[*it].parents.push_back(i);
      
      for (int j=0; j<getNfaceDofs(); j++)
      {
        faces[*it].parentNodes[j][faces[*it].parents.size()-1] = nodeVolumeOrder[getFaceNode(*it,j)];
      }
    }
    nodeVolumeOrder.clear();
  }
  // for each face find its edges
  for (unsigned int i=0; i<faces.size(); i++)
  {
    map<int, int> nodeFaceOrder;
    
    faces[i].children.clear();
    intersection(nodeEdges[faces[i].nodes[0]], nodeEdges[faces[i].nodes[1]], faces[i].children);
    if (create_new && faces[i].children.size() != 1)
    {
      int ne = addEdge(faces[i].nodes[0], faces[i].nodes[1], faces[i].nodes[3]);
      nodeEdges[faces[i].nodes[0]].push_back(ne);
      nodeEdges[faces[i].nodes[1]].push_back(ne);
      nodeEdges[faces[i].nodes[3]].push_back(ne);
      sort(nodeEdges[faces[i].nodes[0]].begin(), nodeEdges[faces[i].nodes[0]].end());
      sort(nodeEdges[faces[i].nodes[1]].begin(), nodeEdges[faces[i].nodes[1]].end());
      sort(nodeEdges[faces[i].nodes[3]].begin(), nodeEdges[faces[i].nodes[3]].end());
      faces[i].children.push_back(ne);
    }
    intersection(nodeEdges[faces[i].nodes[0]], nodeEdges[faces[i].nodes[2]], faces[i].children);
    if (create_new && faces[i].children.size() != 2)
    {
      int ne = addEdge(faces[i].nodes[0], faces[i].nodes[2], faces[i].nodes[5]);
      nodeEdges[faces[i].nodes[0]].push_back(ne);
      nodeEdges[faces[i].nodes[2]].push_back(ne);
      nodeEdges[faces[i].nodes[5]].push_back(ne);
      sort(nodeEdges[faces[i].nodes[0]].begin(), nodeEdges[faces[i].nodes[0]].end());
      sort(nodeEdges[faces[i].nodes[2]].begin(), nodeEdges[faces[i].nodes[2]].end());
      sort(nodeEdges[faces[i].nodes[5]].begin(), nodeEdges[faces[i].nodes[5]].end());
      faces[i].children.push_back(ne);
    }
    intersection(nodeEdges[faces[i].nodes[1]], nodeEdges[faces[i].nodes[2]], faces[i].children);
    if (create_new && faces[i].children.size() != 3)
    {
      int ne = addEdge(faces[i].nodes[1], faces[i].nodes[2], faces[i].nodes[4]);
      nodeEdges[faces[i].nodes[1]].push_back(ne);
      nodeEdges[faces[i].nodes[2]].push_back(ne);
      nodeEdges[faces[i].nodes[4]].push_back(ne);
      sort(nodeEdges[faces[i].nodes[1]].begin(), nodeEdges[faces[i].nodes[1]].end());
      sort(nodeEdges[faces[i].nodes[2]].begin(), nodeEdges[faces[i].nodes[2]].end());
      sort(nodeEdges[faces[i].nodes[4]].begin(), nodeEdges[faces[i].nodes[4]].end());
      faces[i].children.push_back(ne);
    }
    
    for (int j=0; j<getNfaceDofs(); j++) nodeFaceOrder[getFaceNode(i,j)] = j;
    
    for (ChildParentListType::iterator it=faces[i].children.begin(); it!=faces[i].children.end(); it++)
    {
      edges[*it].parents.push_back(i);
      
      for (int j=0; j<getNedgeDofs(); j++)
      {
        edges[*it].parentNodes[j][edges[*it].parents.size()-1] = nodeFaceOrder[getEdgeNode(*it,j)];
      }
    }
    nodeFaceOrder.clear();
  }
  
  nodeFaces.clear();
  nodeEdges.clear();
}





